
library(tidyr)
library(car)
library(robustHD)
library(stargazer)
library(foreign)
library(sandwich)
library(grid)
library(lmtest)
library(ggplot2)
library(plyr)
library(RColorBrewer)
library(cjoint)
library(multiwayvcov)
library(coefplot)

########################################
#### National Survey - Non-Conjoint
########################################

setwd("/Users/Michael/Dropbox/apsr")

#bring in data
socpoc<-read.csv("data/socpocAPSR.csv", stringsAsFactors = F)

#assign ownership groups
renters<-subset(socpoc, own==0)
owners<-subset(socpoc, own==1)

#10% supply support models 
table(renters$supply_dummy)
352/(247+352) #.59
table(owners$supply_dummy)
361/(949+361) #.28


# Table B.2. Support for 10 Percent Supply Increase ####
## bivariate
supply_simple<-(lm(supply_dummy ~ own, socpoc));summary(supply_simple)


# Table B.3. Support for 10 Percent Supply Increase - 7 Point Scale ####
## bivariate
supply_7_simple<-(lm(city_supply ~ own, socpoc));summary(supply_7_simple)
supply_7_simple_se<-sqrt(diag(vcovHC(supply_7_simple, type="HC1")))

# full
supply_7_full<-(lm(city_supply ~ own +scale(ideology)+scale(log(income)) + whitenh  +age + male, subset(socpoc)));summary(supply_7_full)
supply_7_full_se<-sqrt(diag(vcovHC(supply_7_full, type="HC1")))

#full w/ fixed effects
supply_7_full_fe<-(lm(city_supply ~  own +scale(ideology)+ scale(log(income))+ whitenh  + age + male +factor(name), socpoc));summary(supply_7_full_fe)
supply_7_full_fe_se<-sqrt(diag(vcovHC(supply_7_full_fe, type="HC1")))

# Table
stargazer(supply_7_simple,  supply_7_full , supply_7_full_fe, title="Support for 10 Percent Supply Increase - 7 Point Scale", label="supply_7",
          dep.var.labels=c("Support Supply Increase"),dep.var.labels.include = F, dep.var.caption = "",
          column.labels=c("Bivariate","Full","Full with Fixed Effects"),
          covariate.labels=c("Homeownership","Ideology","Income, Log","White, Non-Hispanic","Age","Male"),
          omit.stat = c("ser", "f"), digits=2, align=T, type="text",
          initial.zero = F,  font.size = "small", star.cutoffs = NA, omit.table.layout = "n",
          se=list(supply_7_simple_se, supply_7_full_se, supply_7_full_fe_se), no.space=T,omit=c("name"))


# Figure 6. Renter Support for Supply Citywide, by Average City Rent ####
quantile(socpoc$zri_city, probs=seq(0,1,.2), na.rm=T)
zri_city_values<-c(0,1217,1480,1936,2427,7344)
est1<-rep(NA, length(zri_city_values))
se1<-rep(NA, length(zri_city_values))

for(i in 1:5){
  section<-subset(renters, zri_city>zri_city_values[i] & zri_city<=zri_city_values[i+1])
  est1[i]<-mean(section$supply_dummy,na.rm=T)
  se1[i]<-sqrt((est1[i]*(1-est1[i]))/nrow(section))
}

mod1ests<-as.data.frame(cbind(zri_city_values,est1,se1))
mod1ests$uCI<-est1+se1*1.96
mod1ests$lCI<-est1-se1*1.96

mod1estimates<-mod1ests[1:5,]
colnames(mod1estimates)<-c("Rent","Estimate", "StdErr", "UpperCI", "LowerCI")
mod1estimates$Quintile<-c(1,2,3,4,5)
mod1estimates
mod1estimates$Cost<-factor(c("Least Expensive","2nd","3rd","4th","Most Expensive"))
levels(mod1estimates$Cost)
pd <- position_dodge(0.1)
mod1estimates$Cost <- factor(mod1estimates$Cost, levels = c("Least Expensive","2nd","3rd","4th","Most Expensive"))
mod1estimates
renter_city_supply<-ggplot(mod1estimates, aes(x=Cost, y=Estimate))+scale_y_continuous(limits = c(0, 1))+coord_flip()+
  geom_errorbar(aes(ymin=LowerCI, ymax=UpperCI), width=.1, position=pd) +
  geom_point(position=pd)+  labs(x = "Average City Rent (Quintiles)", y="Support for Supply (%)")+ggtitle("Renters Support for Supply Citywide, by Average Rent") + theme(aspect.ratio = .5)+
  theme_bw()+scale_fill_grey()
pdf("plots/Fig6_renters_city_supply.pdf", width=6, height=3)
print(renter_city_supply)
dev.off()

tiff("plots/Fig6_renters_city_supply.tif", res=600, compression="lzw", height=3, width=6, units="in")
print(renter_city_supply)
dev.off()


# Figure C.8. Renters Suppport for Supply Cityiwde, Average ZIP Rent ####
quantile(socpoc$zri, probs=seq(0,1,.2), na.rm=T)
zri_city_values<-c(0,1204,1526,1958,2488,13000)
est1<-rep(NA, length(zri_city_values))
se1<-rep(NA, length(zri_city_values))

for(i in 1:5){
  section<-subset(renters, zri>zri_city_values[i] & zri<=zri_city_values[i+1])
  est1[i]<-mean(section$supply_dummy,na.rm=T)
  se1[i]<-sqrt((est1[i]*(1-est1[i]))/nrow(section))
}

mod1ests<-as.data.frame(cbind(zri_city_values,est1,se1))
mod1ests$uCI<-est1+se1*1.96
mod1ests$lCI<-est1-se1*1.96

mod1estimates<-mod1ests[1:5,]
colnames(mod1estimates)<-c("Rent","Estimate", "StdErr", "UpperCI", "LowerCI")
mod1estimates$Quintile<-c(1,2,3,4,5)
mod1estimates
mod1estimates$Cost<-factor(c("Least Expensive","2nd","3rd","4th","Most Expensive"))
levels(mod1estimates$Cost)
pd <- position_dodge(0.1)
mod1estimates$Cost <- factor(mod1estimates$Cost, levels = c("Least Expensive","2nd","3rd","4th","Most Expensive"))
mod1estimates
renter_zip_supply<-ggplot(mod1estimates, aes(x=Cost, y=Estimate))+scale_y_continuous(limits = c(0, 1))+coord_flip()+
  geom_errorbar(aes(ymin=LowerCI, ymax=UpperCI), width=.1, position=pd) +
  geom_point(position=pd)+  labs(x = "Average ZIP Rent (Quintiles)", y="Support for Supply (%)")+ggtitle("Renters Support for Supply Citywide, by Average ZIP Rent") + theme(aspect.ratio = .5)+
  theme_bw()+scale_fill_grey()
pdf("plots/FigC8_renters_zip_supply.pdf", width=6, height=3)
print(renter_zip_supply)
dev.off()

tiff("plots/FigC8_renters_zip_supply.tif", res=600, compression="lzw", height=3, width=6, units="in")
print(renter_zip_supply)
dev.off()


# Figure C.9. Homeowners Support for Supply Citywide, Average City Rent ####
quantile(socpoc$zri_city, probs=seq(0,1,.2), na.rm=T)
zri_city_values<-c(0,1217,1480,1936,2427,7344)

est2<-rep(NA, length(zri_city_values))
se2<-rep(NA, length(zri_city_values))

for(i in 1:5){
  section<-subset(owners, zri_city>zri_city_values[i] & zri_city<=zri_city_values[i+1])
  est2[i]<-mean(section$supply_dummy,na.rm=T)
  se2[i]<-sqrt((est2[i]*(1-est2[i]))/nrow(section))
}

mod2ests<-as.data.frame(cbind(zri_city_values,est2,se2))
mod2ests$uCI<-est2+se2*1.96
mod2ests$lCI<-est2-se2*1.96

mod2estimates<-mod2ests[1:5,]
colnames(mod2estimates)<-c("Rent","Estimate", "StdErr", "UpperCI", "LowerCI")
mod2estimates$Quintile<-c(1,2,3,4,5)
mod2estimates
mod2estimates$Cost<-factor(c("Least Expensive","2nd","3rd","4th","Most Expensive"))
levels(mod2estimates$Cost)
pd <- position_dodge(0.1)
mod2estimates$Cost <- factor(mod2estimates$Cost, levels = c("Least Expensive","2nd","3rd","4th","Most Expensive"))
owners_city_supply<-ggplot(mod2estimates, aes(x=Cost, y=Estimate))+scale_y_continuous(limits = c(0, 1))+coord_flip()+
  geom_errorbar(aes(ymin=LowerCI, ymax=UpperCI), width=.1, position=pd) +
  geom_point(position=pd)+ggtitle("Homeowners Support for Supply Citywide, by Average Rent")+
  labs(x = "Average City Rent (Quintiles)", y="Support for Supply (%)")+ theme(aspect.ratio=.5)+
  theme_bw()+scale_fill_grey()
pdf("plots/FigC9_owners_city_supply.pdf", width=6, height=3)
print(owners_city_supply)
dev.off()

tiff("plots/FigC9_owners_city_supply.tif", res=600, compression="lzw", height=3, width=6, units="in")
print(owners_city_supply)
dev.off()


# Supplmentary Figure 3. Homeowners Support for Supply Citywide, Average ZIP Rent ####
quantile(socpoc$zri, probs=seq(0,1,.2), na.rm=T)
zri_city_values<-c(0,1204,1526,1958,2488,13000)

est2<-rep(NA, length(zri_city_values))
se2<-rep(NA, length(zri_city_values))

for(i in 1:5){
  section<-subset(owners, zri>zri_city_values[i] & zri<=zri_city_values[i+1])
  est2[i]<-mean(section$supply_dummy,na.rm=T)
  se2[i]<-sqrt((est2[i]*(1-est2[i]))/nrow(section))
}

mod2ests<-as.data.frame(cbind(zri_city_values,est2,se2))
mod2ests$uCI<-est2+se2*1.96
mod2ests$lCI<-est2-se2*1.96

mod2estimates<-mod2ests[1:5,]
colnames(mod2estimates)<-c("Rent","Estimate", "StdErr", "UpperCI", "LowerCI")
mod2estimates$Quintile<-c(1,2,3,4,5)
mod2estimates
mod2estimates$Cost<-factor(c("Least Expensive","2nd","3rd","4th","Most Expensive"))
levels(mod2estimates$Cost)
pd <- position_dodge(0.1)
mod2estimates$Cost <- factor(mod2estimates$Cost, levels = c("Least Expensive","2nd","3rd","4th","Most Expensive"))
owners_zip_supply<-ggplot(mod2estimates, aes(x=Cost, y=Estimate))+scale_y_continuous(limits = c(0, 1))+coord_flip()+
  geom_errorbar(aes(ymin=LowerCI, ymax=UpperCI), width=.1, position=pd) +
  geom_point(position=pd)+ggtitle("Homeowners Support for Supply Citywide, by Average ZIP Rent")+
  labs(x = "Average ZIP Rent (Quintiles)", y="Support for Supply (%)")+ theme(aspect.ratio=.5)+
  theme_bw()+scale_fill_grey()
pdf("plots/FigSupp3_owners_zip_supply.pdf", width=6, height=3)
print(owners_zip_supply)
dev.off()

tiff("plots/FigSupp3_owners_zip_supply.tif", res=600, compression="lzw", height=3, width=6, units="in")
print(owners_zip_supply)
dev.off()

#### Nimby ban models 
table(renters$ban_dummy)
234/(234+434) #.35
table(owners$ban_dummy)
593/(811+593) #.42


# Table B.4. Support for Neighborhood Ban ####
## bivariate
ban_simple<-(lm(ban_dummy ~ own, socpoc));summary(ban_simple)
ban_simple_se<-sqrt(diag(vcovHC(ban_simple, type="HC1")))

# full
ban_full<-(lm(ban_dummy ~ own +scale(ideology)+scale(log(income)) + whitenh  +age + male, socpoc));summary(ban_full)
ban_full_se<-sqrt(diag(vcovHC(ban_full, type="HC1")))

#full w/ fixed effects
ban_full_fe<-(lm(ban_dummy ~  own +scale(ideology)+ scale(log(income))+ whitenh  + age + male +factor(name), socpoc));summary(ban_full_fe)
ban_full_fe_se<-sqrt(diag(vcovHC(ban_full_fe, type="HC1")))

# Table
stargazer(ban_simple,  ban_full , ban_full_fe, title="Support for Ban on Neighborhood Development", label="ban_dummy",
          dep.var.labels=c("Support NIMBY Ban"),dep.var.labels.include = F, dep.var.caption = "",
          column.labels=c("Bivariate","Full","Full with Fixed Effects"),
          covariate.labels=c("Homeownership","Ideology","Income, Log","White, Non-Hispanic","Age","Male"),
          omit.stat = c("ser", "f"), digits=2, align=T, type="text",
          initial.zero = F,  font.size = "small",star.cutoffs = NA, omit.table.layout = "n",
          se=list(ban_simple_se, ban_full_se, ban_full_fe_se), no.space=T, omit=c("name"))

# Table B.5. Support for Neighborhood Ban 7 point scale ####
#simple
ban_simple<-(lm(neighborhood_ban ~ own, socpoc));summary(ban_simple)
ban_simple_se<-sqrt(diag(vcovHC(ban_simple, type="HC1")))

# full
ban_full<-(lm(neighborhood_ban ~ own +scale(ideology)+scale(log(income)) + whitenh  +age + male, socpoc));summary(ban_full)
ban_full_se<-sqrt(diag(vcovHC(ban_full, type="HC1")))

#full w/ fixed effects
ban_full_fe<-(lm(neighborhood_ban ~  own +scale(ideology)+ scale(log(income))+ whitenh  + age + male +factor(name), socpoc));summary(ban_full_fe)
ban_full_fe_se<-sqrt(diag(vcovHC(ban_full_fe, type="HC1")))

# Table
stargazer(ban_simple,  ban_full , ban_full_fe, title="Support for Ban on Neighborhood Development - 7 Point Scale", label="neighborhood_ban_7",
          dep.var.labels=c("Support NIMBY Ban"),dep.var.labels.include = F, dep.var.caption = "",
          column.labels=c("Bivariate","Full","Full with Fixed Effects"),
          covariate.labels=c("Homeownership","Ideology","Income, Log","White, Non-Hispanic","Age","Male"),
          omit.stat = c("ser", "f"), digits=2, align=T, type="text",
          initial.zero = F,  font.size = "small",star.cutoffs = NA, omit.table.layout = "n",
          se=list(ban_simple_se, ban_full_se, ban_full_fe_se), no.space=T, omit=c("name"))


####################################
#### National Survey - Conjoint 
####################################

setwd("/Users/Michael/Dropbox/apsr")

conjoint4<-read.csv("data/conjointDataAPSR.csv")
table(conjoint4$own)
#relevel
conjoint4$distance <- factor(conjoint4$distance,levels= c("2 miles (40 minute walk)", "1 mile (20 minute walk)", "1/2 mile (10 minute walk)","1/8 mile (2 minute walk)"))
conjoint4$community <- factor(conjoint4$community,levels= c("No opinion", "Support the building", "Oppose the building"))
conjoint4$affordable <- factor(conjoint4$affordable,levels= c("None of the units", "One-quarter of the units", "Half of the units", "All of the units"))
conjoint4$height <- factor(conjoint4$height,levels= c("2 stories", "3 stories", "6 stories", "12 stories"))
conjoint4$site <- factor(conjoint4$site, levels=c("Empty building","Parking lot","Open field","Historically-designated building"))
names(conjoint4)
#reclassify items as factor
cols<-c("own", "whitenh", "nearby","conjoint_first","rich", "luxury")

conjoint4[cols] <- data.frame(apply(conjoint4[c(cols)], 2, as.factor))

# dummies
conjoint4$liberal<-as.factor(ifelse(conjoint4$ideology>4,1,ifelse(conjoint4$ideology<4,0, NA)))
conjoint4$city_interest_low<-as.factor(ifelse(conjoint4$city_interest<0,1,0))

#define subgroups/dummies
renters<-subset(conjoint4, own==0)
owners<-subset(conjoint4, own==1)

#define affordability levels
renters_aff<-subset(renters, aff_housing==1)
renters_lux<-subset(renters, aff_housing==0)

owners_aff<-subset(owners, aff_housing==1)
owners_lux<-subset(owners, aff_housing==0)


# FIGURE 3. Homeowners, Proximity by Affordability ####
owners_luxury_mod<-amce(select ~ distance  + community + height + site + tenant + units,
                        data= owners_lux , cluster=T,  respondent.id = "CaseID")
owners_affordable_mod<-amce(select ~ distance  + community + height + site + tenant + units,
                            data=owners_aff, cluster=T,  respondent.id = "CaseID")
owners_luxury_mod_frame<-data.frame(Variable=(summary(owners_luxury_mod)$amce)$Level,
                                    Coefficient = (summary(owners_luxury_mod)$amce)$Estimate,
                                    SE=(summary(owners_luxury_mod)$amce)$'Std. Err',
                                    modelName="Market Rate")
owners_affordable_mod_frame<-data.frame(Variable=(summary(owners_affordable_mod)$amce)$Level,
                                        Coefficient = (summary(owners_affordable_mod)$amce)$Estimate,
                                        SE=(summary(owners_affordable_mod)$amce)$'Std. Err',
                                        modelName="Affordable")
ownersPriceFrame<-data.frame(rbind(owners_luxury_mod_frame, owners_affordable_mod_frame))
ownersPriceFrame<-subset(ownersPriceFrame, Variable=="1/8 mile (2 minute walk)"|Variable=="1/2 mile (10 minute walk)"|Variable=="1 mile (20 minute walk)")
ownersPriceFrameIntercepts<-data.frame(Variable=c("2 miles (40 minute walk)", "2 miles (40 minute walk)") ,Coefficient=c(0,0), SE=c(0,0), modelName=c("Market Rate", "Affordable"))
ownersPriceFrame<-data.frame(rbind(ownersPriceFrame,ownersPriceFrameIntercepts))
interval1<-qnorm((1-.9)/2)
interval2<-qnorm((1-.95)/2)
ownersPriceFrame$Variable <- factor(ownersPriceFrame$Variable, levels = c("1/8 mile (2 minute walk)","1/2 mile (10 minute walk)", "1 mile (20 minute walk)", "2 miles (40 minute walk)"))
owners_price_nimby<-ggplot(ownersPriceFrame, aes(colour=modelName, shape=modelName))+ scale_y_continuous(limits = c(-.20, .20))
owners_price_nimby<-owners_price_nimby+theme_bw()+scale_colour_grey(end=.5)+geom_hline(yintercept=0, colour=gray(1/2), lty=2)
owners_price_nimby<-owners_price_nimby+geom_linerange(aes(x=Variable, ymin=Coefficient-SE*interval1, 
                                                          ymax=Coefficient+SE*interval1), lwd=1, position=position_dodge(width=1/2))
owners_price_nimby<-owners_price_nimby+geom_pointrange(aes(x=Variable, y=Coefficient, ymin=Coefficient-SE*interval2,
                                                           ymax=Coefficient+SE*interval2), lwd=1/2,
                                                       position=position_dodge(width=1/2), fill="WHITE")
owners_price_nimby<-owners_price_nimby+coord_flip()+labs(y="Change in Probability Building Preferred")
owners_price_nimby<-owners_price_nimby+theme(legend.title=element_blank(),axis.title.y=element_blank())+ theme(aspect.ratio = .5)
owners_price_nimby<-owners_price_nimby+theme(plot.margin=unit(c(0,0,0,0),"mm"))+ggtitle("Homeowners, Proximity by Affordability") + guides(colour=guide_legend(reverse=T), shape=guide_legend(reverse=T))
pdf("plots/Fig3_owners_price_nimby.pdf", width=6, height=3)
print(owners_price_nimby)
dev.off()

tiff("plots/Fig3_owners_price_nimby.tif", res=600, compression="lzw", height=3, width=6, units="in")
print(owners_price_nimby)
dev.off()


# FIGURE 4. Renters, Proximity by Affordability ####
renters_luxury_mod<-amce(select ~ distance  + community + height + site + tenant + units,
                         data= renters_lux , cluster=T,  respondent.id = "CaseID")
renters_affordable_mod<-amce(select ~ distance  + community + height + site + tenant + units,
                             data=renters_aff, cluster=T,  respondent.id = "CaseID")
renters_luxury_mod_frame<-data.frame(Variable=(summary(renters_luxury_mod)$amce)$Level,
                                     Coefficient = (summary(renters_luxury_mod)$amce)$Estimate,
                                     SE=(summary(renters_luxury_mod)$amce)$'Std. Err',
                                     modelName="Market Rate")
renters_affordable_mod_frame<-data.frame(Variable=(summary(renters_affordable_mod)$amce)$Level,
                                         Coefficient = (summary(renters_affordable_mod)$amce)$Estimate,
                                         SE=(summary(renters_affordable_mod)$amce)$'Std. Err',
                                         modelName="Affordable")
rentersPriceFrame<-data.frame(rbind(renters_luxury_mod_frame, renters_affordable_mod_frame))
rentersPriceFrame<-subset(rentersPriceFrame, Variable=="1/8 mile (2 minute walk)"|Variable=="1/2 mile (10 minute walk)"|Variable=="1 mile (20 minute walk)")
rentersPriceFrameIntercepts<-data.frame(Variable=c("2 miles (40 minute walk)", "2 miles (40 minute walk)") ,Coefficient=c(0,0), SE=c(0,0), modelName=c("Market Rate", "Affordable"))
rentersPriceFrame<-data.frame(rbind(rentersPriceFrame,rentersPriceFrameIntercepts))
interval1<-qnorm((1-.9)/2)
interval2<-qnorm((1-.95)/2)
rentersPriceFrame$Variable <- factor(rentersPriceFrame$Variable, levels = c("1/8 mile (2 minute walk)","1/2 mile (10 minute walk)", "1 mile (20 minute walk)", "2 miles (40 minute walk)"))
renters_price_nimby<-ggplot(rentersPriceFrame, aes(colour=modelName, shape=modelName)) + scale_y_continuous(limits = c(-.20, .20))
renters_price_nimby<-renters_price_nimby+theme_bw()+scale_colour_grey(end=.5)+geom_hline(yintercept=0, colour=gray(1/2), lty=2)
renters_price_nimby<-renters_price_nimby+geom_linerange(aes(x=Variable, ymin=Coefficient-SE*interval1, 
                                                            ymax=Coefficient+SE*interval1), lwd=1, position=position_dodge(width=1/2))
renters_price_nimby<-renters_price_nimby+geom_pointrange(aes(x=Variable, y=Coefficient, ymin=Coefficient-SE*interval2,
                                                             ymax=Coefficient+SE*interval2), lwd=1/2,
                                                         position=position_dodge(width=1/2), fill="WHITE")
renters_price_nimby<-renters_price_nimby+coord_flip()+ labs(y="Change in Probability Building Preferred")
renters_price_nimby<-renters_price_nimby+theme(legend.title=element_blank(), axis.title.y=element_blank())+ theme(aspect.ratio = .5)
renters_price_nimby<-renters_price_nimby+theme(plot.margin=unit(c(0,0,0,0),"mm")) +ggtitle("Renters, Proximity by Affordability")+ guides(colour=guide_legend(reverse=T), shape=guide_legend(reverse=T))
pdf("plots/Fig4_renters_price_nimby.pdf", width=6, height=3)
print(renters_price_nimby)
dev.off()

tiff("plots/Fig4_renters_price_nimby.tif", res=600, compression="lzw", height=3, width=6, units="in")
print(renters_price_nimby)
dev.off()


# FIGURE 5, Renters, Nimby by Affordability, Quintile City####
quantile(conjoint4$zri_city, probs=seq(0,1,.1), na.rm=T) # define quintiles

zri_city_values<-c(0,1217,1480,1936,2427,7344)
est1<-rep(NA, length(zri_city_values))
se1<-rep(NA, length(zri_city_values))
for(i in 1:5){
  mod1<-amce(select ~ distance+ community + height + site + tenant + units,
             data=subset(renters, zri_city>zri_city_values[i] & zri_city<=zri_city_values[i+1]&luxury==1), cluster=T, respondent.id = "CaseID")
  est1[i]<-summary(mod1)$amce[5,3]
  se1[i]<-summary(mod1)$amce[5,4]
}
mod1ests<-as.data.frame(cbind(zri_city_values,est1,se1))
mod1ests$uCI<-est1+se1*1.96
mod1ests$lCI<-est1-se1*1.96
est2<-rep(NA, length(zri_city_values))
se2<-rep(NA, length(zri_city_values))
for(i in 1:5){
  mod2<-amce(select ~ distance+ community + height + site + tenant + units,
             data=subset(renters, zri_city>zri_city_values[i] & zri_city<=zri_city_values[i+1]&luxury==0), cluster=T, respondent.id = "CaseID")
  est2[i]<-summary(mod2)$amce[5,3]
  se2[i]<-summary(mod2)$amce[5,4]
}
mod2ests<-as.data.frame(cbind(zri_city_values,est2,se2))
mod2ests$uCI<-est2+se2*1.96
mod2ests$lCI<-est2-se2*1.96
#combine data
mod1<-mod1ests[-6,]
mod1$quintle<-c("Least Expensive","2nd","3rd","4th","Most Expensive")
mod1$modelName<-"Market Rate"
names(mod1)<-c("zri","est","se","uci","lci","Quintile","modelName")
mod2<-mod2ests[-6,]
mod2$quintle<-c("Least Expensive","2nd","3rd","4th","Most Expensive")
mod2$modelName<-"Affordable"
names(mod2)<-c("zri","est","se","uci","lci","Quintile","modelName")
modFrame<-data.frame(rbind(mod1,mod2))
modFrame$Quintile <- factor(modFrame$Quintile, levels = c("Least Expensive", "2nd","3rd", "4th", "Most Expensive"))
interval1<-qnorm((1-.9)/2)
interval2<-qnorm((1-.95)/2)
modFrame
modFrame$modelName<-factor(modFrame$modelName, levels=c("Market Rate","Affordable"))
renters_type_nimby<-ggplot(modFrame, aes(colour=modelName, shape=modelName))+ scale_y_continuous(limits = c(-.25, .25))
renters_type_nimby<-renters_type_nimby+theme_bw()+scale_colour_grey(end=.5)+geom_hline(yintercept=0, colour=gray(1/2), lty=2)
renters_type_nimby<-renters_type_nimby+geom_linerange(aes(x=Quintile, ymin=est-se*interval1, 
                                                          ymax=est+se*interval1), lwd=1, position=position_dodge(width=1/2))
#+scale_color_manual(values=c('#F8766D', '#00BFC4'))

renters_type_nimby<-renters_type_nimby+geom_pointrange(aes(x=Quintile, y=est, ymin=est-se*interval2,
                                                           ymax=est+se*interval2), lwd=1/2,
                                                       position=position_dodge(width=1/2),  fill="WHITE")
renters_type_nimby<-renters_type_nimby+coord_flip()+ labs(y="Change in Probability Building Preferred",x="Average Rent by City (Quintiles)")
renters_type_nimby<-renters_type_nimby+theme(legend.title=element_blank()) + theme(aspect.ratio = .5)
renters_type_nimby<-renters_type_nimby+theme(plot.margin=unit(c(0,0,0,0),"mm"))+ ggtitle("Renters, Proximity by Affordability and Average Rent by City") +guides(colour=guide_legend(reverse=T), shape=guide_legend(reverse=T))
pdf("plots/Fig5_renters_type_nimby.pdf", width=6, height=3)
print(renters_type_nimby)
dev.off()

tiff("plots/Fig5_renters_type_nimby.tif", res=600, compression="lzw", height=3, width=6, units="in")
print(renters_type_nimby)
dev.off()

# FIGURE 7, Renters, nimby by price anxiety ####
renters_city_low_mod<-amce(select ~ distance  + community + height + site + tenant + units,
                           data=subset(renters_lux, city_interest<0) , cluster=T,  respondent.id = "CaseID")
renters_city_high_mod<-amce(select ~ distance  + community + height + site + tenant + units,
                            data=subset(renters_lux, city_interest>=0), cluster=T,  respondent.id = "CaseID")
renters_city_low_mod_frame<-data.frame(Variable=(summary(renters_city_low_mod)$amce)$Level,
                                       Coefficient = (summary(renters_city_low_mod)$amce)$Estimate,
                                       SE=(summary(renters_city_low_mod)$amce)$'Std. Err',
                                       modelName="Price Anxious")
renters_city_high_mod_frame<-data.frame(Variable=(summary(renters_city_high_mod)$amce)$Level,
                                        Coefficient = (summary(renters_city_high_mod)$amce)$Estimate,
                                        SE=(summary(renters_city_high_mod)$amce)$'Std. Err',
                                        modelName="Price Neutral")
rentersCityFrame<-data.frame(rbind(renters_city_low_mod_frame, renters_city_high_mod_frame))
rentersCityFrame<-subset(rentersCityFrame, Variable=="1/8 mile (2 minute walk)"|Variable=="1/2 mile (10 minute walk)"|Variable=="1 mile (20 minute walk)")
rentersCityFrameIntercepts<-data.frame(Variable=c("2 miles (40 minute walk)", "2 miles (40 minute walk)") ,Coefficient=c(0,0), SE=c(0,0), modelName=c("Price Anxious", "Price Neutral"))
rentersCityFrame<-data.frame(rbind(rentersCityFrame,rentersCityFrameIntercepts))
interval1<-qnorm((1-.9)/2)
interval2<-qnorm((1-.95)/2)
rentersCityFrame$Variable <- factor(rentersCityFrame$Variable, levels = c("1/8 mile (2 minute walk)","1/2 mile (10 minute walk)", "1 mile (20 minute walk)", "2 miles (40 minute walk)"))
renters_anxious_nimby<-ggplot(rentersCityFrame, aes(colour=modelName, shape=modelName))+ scale_y_continuous(limits = c(-.25, .25))
renters_anxious_nimby<-renters_anxious_nimby+theme_bw()+scale_colour_grey(end=.5)+geom_hline(yintercept=0, colour=gray(1/2), lty=2)+scale_shape_manual(values=c(4,5))
#+ scale_color_manual(values=c('#990099','#33CC00'))
renters_anxious_nimby<-renters_anxious_nimby+geom_linerange(aes(x=Variable, ymin=Coefficient-SE*interval1, 
                                                                ymax=Coefficient+SE*interval1), lwd=1, position=position_dodge(width=1/2))
renters_anxious_nimby<-renters_anxious_nimby+geom_pointrange(aes(x=Variable, y=Coefficient, ymin=Coefficient-SE*interval2,
                                                                 ymax=Coefficient+SE*interval2), lwd=1/2,
                                                             position=position_dodge(width=1/2), fill="WHITE")
renters_anxious_nimby<-renters_anxious_nimby+coord_flip()+ labs(y="Change in Probability Building Preferred")
renters_anxious_nimby<-renters_anxious_nimby+theme(legend.title=element_blank(),axis.title.y=element_blank())+theme(aspect.ratio = .5)
renters_anxious_nimby<-renters_anxious_nimby+theme(plot.margin=unit(c(0,0,0,0),"mm")) + ggtitle("Renters, Proximity by Price Anxiety (Market Rate)")+ theme(aspect.ratio = .5)+ guides(colour=guide_legend(reverse=T), shape=guide_legend(reverse=T))
pdf("plots/Fig7_renters_anxious_nimby.pdf", width=6, height=3)
print(renters_anxious_nimby)
dev.off()

tiff("plots/Fig7_renters_anxious_nimby.tif", res=600, compression="lzw", height=3, width=6, units="in")
print(renters_anxious_nimby)
dev.off()


# FIGURE C.7. Renters, nimby by price anxiety to affordable ####
renters_city_low_mod<-amce(select ~ distance  + community + height + site + tenant + units,
                           data=subset(renters_aff, city_interest<0) , cluster=T,  respondent.id = "CaseID")
renters_city_high_mod<-amce(select ~ distance  + community + height + site + tenant + units,
                            data=subset(renters_aff, city_interest>=0), cluster=T,  respondent.id = "CaseID")
renters_city_low_mod_frame<-data.frame(Variable=(summary(renters_city_low_mod)$amce)$Level,
                                       Coefficient = (summary(renters_city_low_mod)$amce)$Estimate,
                                       SE=(summary(renters_city_low_mod)$amce)$'Std. Err',
                                       modelName="Price Anxious")
renters_city_high_mod_frame<-data.frame(Variable=(summary(renters_city_high_mod)$amce)$Level,
                                        Coefficient = (summary(renters_city_high_mod)$amce)$Estimate,
                                        SE=(summary(renters_city_high_mod)$amce)$'Std. Err',
                                        modelName="Price Neutral")
rentersCityFrame<-data.frame(rbind(renters_city_low_mod_frame, renters_city_high_mod_frame))
rentersCityFrame<-subset(rentersCityFrame, Variable=="1/8 mile (2 minute walk)"|Variable=="1/2 mile (10 minute walk)"|Variable=="1 mile (20 minute walk)")
rentersCityFrameIntercepts<-data.frame(Variable=c("2 miles (40 minute walk)", "2 miles (40 minute walk)") ,Coefficient=c(0,0), SE=c(0,0), modelName=c("Price Anxious", "Price Neutral"))
rentersCityFrame<-data.frame(rbind(rentersCityFrame,rentersCityFrameIntercepts))
interval1<-qnorm((1-.9)/2)
interval2<-qnorm((1-.95)/2)
rentersCityFrame$Variable <- factor(rentersCityFrame$Variable, levels = c("1/8 mile (2 minute walk)","1/2 mile (10 minute walk)", "1 mile (20 minute walk)", "2 miles (40 minute walk)"))
renters_anxious_nimby<-ggplot(rentersCityFrame, aes(colour=modelName, shape=modelName))+ scale_y_continuous(limits = c(-.25, .25))
renters_anxious_nimby<-renters_anxious_nimby+theme_bw()+scale_colour_grey(end=.5)+geom_hline(yintercept=0, colour=gray(1/2), lty=2)+scale_shape_manual(values=c(4,5))
#+  scale_color_manual(values=c('#990099','#33CC00'))
renters_anxious_nimby<-renters_anxious_nimby+geom_linerange(aes(x=Variable, ymin=Coefficient-SE*interval1, 
                                                                ymax=Coefficient+SE*interval1), lwd=1, position=position_dodge(width=1/2))
renters_anxious_nimby<-renters_anxious_nimby+geom_pointrange(aes(x=Variable, y=Coefficient, ymin=Coefficient-SE*interval2,
                                                                 ymax=Coefficient+SE*interval2), lwd=1/2,
                                                             position=position_dodge(width=1/2), fill="WHITE")
renters_anxious_nimby<-renters_anxious_nimby+coord_flip()+ labs(y="Change in Probability Building Preferred")
renters_anxious_nimby<-renters_anxious_nimby+theme(legend.title=element_blank(),axis.title.y=element_blank())+theme(aspect.ratio = .5)
renters_anxious_nimby<-renters_anxious_nimby+theme(plot.margin=unit(c(0,0,0,0),"mm")) + ggtitle("Renters, Proximity by Price Anxiety (Affordable)")+ theme(aspect.ratio = .5)+ guides(colour=guide_legend(reverse=T), shape=guide_legend(reverse=T))
pdf("plots/FigC7_renters_anxious_nimby_affordable.pdf", width=6, height=3)
print(renters_anxious_nimby)
dev.off()

tiff("plots/FigC7_renters_anxious_nimby_affordable.tif", res=600, compression="lzw", height=3, width=6, units="in")
print(renters_anxious_nimby)
dev.off()


# Supplymentary Index, Figure 2. Renters rent burden, employed, citywide ####
renters_employed<-subset(renters, employed==1)
quantile(renters_employed$burden_city, probs=seq(0,1,.2), na.rm=T)
burden_values<-c(.05,.33,.48,.68,1.03,12)
est1<-rep(NA, length(burden_values))
se1<-rep(NA, length(burden_values))
for(i in 1:5){
  mod1<-amce(select ~ distance+ community + height + site + tenant + units,
             data=subset(renters_employed, burden_city>burden_values[i] & burden_city<=burden_values[i+1]&luxury==1), cluster=T, respondent.id = "CaseID")
  est1[i]<-summary(mod1)$amce[5,3]
  se1[i]<-summary(mod1)$amce[5,4]
}
mod1ests<-as.data.frame(cbind(burden_values,est1,se1))
mod1ests$uCI<-est1+se1*1.96
mod1ests$lCI<-est1-se1*1.96
est2<-rep(NA, length(burden_values))
se2<-rep(NA, length(burden_values))
for(i in 1:5){
  mod2<-amce(select ~ distance+ community + height + site + tenant + units,
             data=subset(renters_employed, burden_city>burden_values[i] & burden_city<=burden_values[i+1]&luxury==0), cluster=T, respondent.id = "CaseID")
  est2[i]<-summary(mod2)$amce[5,3]
  se2[i]<-summary(mod2)$amce[5,4]
}
mod2ests<-as.data.frame(cbind(burden_values,est2,se2))
mod2ests$uCI<-est2+se2*1.96
mod2ests$lCI<-est2-se2*1.96
#combine data
mod1<-mod1ests[-6,]
mod1$quintle<-c("Least Burdened","2nd","3rd","4th","Most Burdened")
mod1$modelName<-"Market Rate"
names(mod1)<-c("zri","est","se","uci","lci","Quintile","modelName")
mod2<-mod2ests[-6,]
mod2$quintle<-c("Least Burdened","2nd","3rd","4th","Most Burdened")
mod2$modelName<-"Affordable"
names(mod2)<-c("zri","est","se","uci","lci","Quintile","modelName")
modFrame<-data.frame(rbind(mod1,mod2))
modFrame$Quintile <- factor(modFrame$Quintile, levels = c("Least Burdened", "2nd","3rd", "4th", "Most Burdened"))
interval1<-qnorm((1-.9)/2)
interval2<-qnorm((1-.95)/2)
modFrame
modFrame$modelName<-factor(modFrame$modelName, levels=c("Market Rate","Affordable"))
renters_type_nimby<-ggplot(modFrame, aes(colour=modelName, shape=modelName))+ scale_y_continuous(limits = c(-.5, .5))
renters_type_nimby<-renters_type_nimby+theme_bw()+scale_colour_grey(end=.5)+geom_hline(yintercept=0, colour=gray(1/2), lty=2)
renters_type_nimby<-renters_type_nimby+geom_linerange(aes(x=Quintile, ymin=est-se*interval1, 
                                                          ymax=est+se*interval1), lwd=1, position=position_dodge(width=1/2))
renters_type_nimby<-renters_type_nimby+geom_pointrange(aes(x=Quintile, y=est, ymin=est-se*interval2,
                                                           ymax=est+se*interval2), lwd=1/2,
                                                       position=position_dodge(width=1/2),  fill="WHITE")
renters_type_nimby<-renters_type_nimby+coord_flip()+ labs(y="Change in Probability Building Preferred",x="City Average Rent, by Quintile")
renters_type_nimby<-renters_type_nimby+theme(legend.title=element_blank())+ theme(aspect.ratio = .5)
renters_type_nimby<-renters_type_nimby+theme(plot.margin=unit(c(0,0,0,0),"mm"))+ggtitle("Renters, Proximity by Affordability and Rent Burden by City")+ guides(colour=guide_legend(reverse=T), shape=guide_legend(reverse=T))
pdf("plots/FigSupp2_renters_quintile_burden_city.pdf", width=6, height=3)
print(renters_type_nimby)
dev.off()

tiff("plots/FigSupp2_renters_quintile_burden_city.tif", res=600, compression="lzw", height=3, width=6, units="in")
print(renters_type_nimby)
dev.off()


# Figure C.1. Homeowner Proximity by Income ####
summary(owners$income) #median==80,000
owners_liberal<-amce(select ~ distance + affordable + community + height + site + tenant + units,
                     data=subset(owners, income>80000), cluster=T,  respondent.id = "CaseID")
owners_conservatives<-amce(select ~ distance + affordable + community + height + site + tenant + units,
                           data=subset(owners, income<=80000), cluster=T,  respondent.id = "CaseID")
owners_liberal_frame<-data.frame(Variable=(summary(owners_liberal)$amce)$Level,
                                 Coefficient = (summary(owners_liberal)$amce)$Estimate,
                                 SE=(summary(owners_liberal)$amce)$'Std. Err',
                                 modelName="Below Median Income")
owners_liberal_frame
owners_conservatives_frame<-data.frame(Variable=(summary(owners_conservatives)$amce)$Level,
                                       Coefficient = (summary(owners_conservatives)$amce)$Estimate,
                                       SE=(summary(owners_conservatives)$amce)$'Std. Err',
                                       modelName="Above Median Income")
ideologyFrame<-data.frame(rbind(owners_liberal_frame, owners_conservatives_frame))
ideologyFrame<-subset(ideologyFrame, Variable=="1/8 mile (2 minute walk)"|Variable=="1/2 mile (10 minute walk)"|Variable=="1 mile (20 minute walk)")
ideologyFrameIntercepts<-data.frame(Variable=c("2 miles (40 minute walk)", "2 miles (40 minute walk)") ,Coefficient=c(0,0), SE=c(0,0), modelName=c("Below Median Income", "Above Median Income"))
ideologyFrame<-data.frame(rbind(ideologyFrame,ideologyFrameIntercepts))
interval1<-qnorm((1-.9)/2)
interval2<-qnorm((1-.95)/2)
ideologyFrame$Variable <- factor(ideologyFrame$Variable, levels = c("1/8 mile (2 minute walk)","1/2 mile (10 minute walk)", "1 mile (20 minute walk)", "2 miles (40 minute walk)"))
ideology_affordable<-ggplot(ideologyFrame, aes(colour=modelName, shape=modelName))+ scale_y_continuous(limits = c(-.25, .25))
ideology_affordable<-ideology_affordable+theme_bw()+scale_colour_grey(end=.5)+geom_hline(yintercept=0, colour=gray(1/2), lty=2)
ideology_affordable<-ideology_affordable+geom_linerange(aes(x=Variable, ymin=Coefficient-SE*interval1, 
                                                            ymax=Coefficient+SE*interval1), lwd=1, position=position_dodge(width=1/2))
ideology_affordable<-ideology_affordable+geom_pointrange(aes(x=Variable, y=Coefficient, ymin=Coefficient-SE*interval2,
                                                             ymax=Coefficient+SE*interval2), lwd=1/2,
                                                         position=position_dodge(width=1/2),  fill="WHITE")
ideology_affordable<-ideology_affordable+coord_flip()+ labs(y="Change in Probability Building Preferred")
ideology_affordable<-ideology_affordable+theme(legend.title=element_blank(),axis.title.y=element_blank())+ theme(aspect.ratio = .5)
ideology_affordable<-ideology_affordable+theme(plot.margin=unit(c(0,0,0,0),"mm"))+ggtitle("Homeowners, Proximity by Income")+guides(colour=guide_legend(reverse=T), shape=guide_legend(reverse=T))
pdf("plots/FigC1_owners_nimby_income.pdf", width=6, height=3)
print(ideology_affordable)
dev.off()

tiff("plots/FigC1_owners_nimby_income.tif", res=600, compression="lzw", height=3, width=6, units="in")
print(ideology_affordable)
dev.off()

# Figure C.2 Homeowner Proximity by Ideology ####
owners_liberal<-amce(select ~ distance + affordable  + community + height + site + tenant + units,
                     data=subset(owners_aff, ideology>4), cluster=T,  respondent.id = "CaseID")
owners_conservatives<-amce(select ~ distance + affordable  + community + height + site + tenant + units,
                           data=subset(owners_aff, ideology<4), cluster=T,  respondent.id = "CaseID")
owners_liberal_frame<-data.frame(Variable=(summary(owners_liberal)$amce)$Level,
                                 Coefficient = (summary(owners_liberal)$amce)$Estimate,
                                 SE=(summary(owners_liberal)$amce)$'Std. Err',
                                 modelName="Liberals")
owners_liberal_frame
owners_conservatives_frame<-data.frame(Variable=(summary(owners_conservatives)$amce)$Level,
                                       Coefficient = (summary(owners_conservatives)$amce)$Estimate,
                                       SE=(summary(owners_conservatives)$amce)$'Std. Err',
                                       modelName="Conservatives")
ideologyFrame<-data.frame(rbind( owners_conservatives_frame,owners_liberal_frame))
ideologyFrame<-subset(ideologyFrame, Variable=="1/8 mile (2 minute walk)"|Variable=="1/2 mile (10 minute walk)"|Variable=="1 mile (20 minute walk)")
ideologyFrameIntercepts<-data.frame(Variable=c("2 miles (40 minute walk)", "2 miles (40 minute walk)") ,Coefficient=c(0,0), SE=c(0,0), modelName=c( "Conservatives", "Liberals"))
ideologyFrame<-data.frame(rbind(ideologyFrame,ideologyFrameIntercepts))
interval1<-qnorm((1-.9)/2)
interval2<-qnorm((1-.95)/2)
ideologyFrame$Variable <- factor(ideologyFrame$Variable, levels = c("1/8 mile (2 minute walk)","1/2 mile (10 minute walk)", "1 mile (20 minute walk)", "2 miles (40 minute walk)"))
ideology_affordable<-ggplot(ideologyFrame, aes(colour=modelName, shape=modelName))+ scale_y_continuous(limits = c(-.25, .25))
ideology_affordable<-ideology_affordable+theme_bw()+scale_colour_grey(end=.5)+geom_hline(yintercept=0, colour=gray(1/2), lty=2)
ideology_affordable<-ideology_affordable+geom_linerange(aes(x=Variable, ymin=Coefficient-SE*interval1, 
                                                            ymax=Coefficient+SE*interval1), lwd=1, position=position_dodge(width=1/2))
ideology_affordable<-ideology_affordable+geom_pointrange(aes(x=Variable, y=Coefficient, ymin=Coefficient-SE*interval2,
                                                             ymax=Coefficient+SE*interval2), lwd=1/2,
                                                         position=position_dodge(width=1/2),  fill="WHITE")
ideology_affordable<-ideology_affordable+coord_flip()+ labs(y="Change in Probability Building Preferred")
ideology_affordable<-ideology_affordable+theme(legend.title=element_blank(),axis.title.y=element_blank())+ theme(aspect.ratio = .5)
ideology_affordable<-ideology_affordable+theme(plot.margin=unit(c(0,0,0,0),"mm"))+ggtitle("Homeowners, Proximity by Ideology")+guides(colour=guide_legend(reverse=T), shape=guide_legend(reverse=T))
pdf("plots/FigC2_owners_nimby_ideology.pdf", width=6, height=3)
print(ideology_affordable)
dev.off()

tiff("plots/FigC2_owners_nimby_ideology.tif", res=600, compression="lzw", height=3, width=6, units="in")
print(ideology_affordable)
dev.off()

# Figure C.3 Proximity by Affordability Homeowners and Renters ####
affordable_values<-c("None of the units", "One-quarter of the units", "Half of the units", "All of the units")
est1<-rep(NA, length(affordable_values))
se1<-rep(NA, length(affordable_values))
for(i in 1:4){
  mod1<-amce(select ~ distance + community + height + site + tenant + units,
             data=subset(owners, affordable==affordable_values[i]), cluster=T, respondent.id = "CaseID")
  est1[i]<-summary(mod1)$amce[5,3]
  se1[i]<-summary(mod1)$amce[5,4]
}
mod1ests<-as.data.frame(cbind(affordable_values,est1,se1))
mod1ests$uCI<-est1+se1*1.96
mod1ests$lCI<-est1-se1*1.96

est2<-rep(NA, length(affordable_values))
se2<-rep(NA, length(affordable_values))
for(i in 1:4){
  mod2<-amce(select ~ distance + community + height + site + tenant + units,
             data=subset(renters, affordable==affordable_values[i]), cluster=T, respondent.id = "CaseID")
  est2[i]<-summary(mod2)$amce[5,3]
  se2[i]<-summary(mod2)$amce[5,4]
}
mod2ests<-as.data.frame(cbind(affordable_values,est2,se2))
mod2ests$uCI<-est1+se1*1.96
mod2ests$lCI<-est1-se1*1.96
mod2ests
#combine data

mod1<-mod1ests
class(mod1$est)
mod1$quintle<-c("None of the units", "One-quarter of the units", "Half of the units", "All of the units")
mod1$modelName<-"Homeowners"
names(mod1)<-c("affordability","est","se","uci","lci","Quintile","modelName")
mod2<-mod2ests
mod2$quintle<-c("None of the units", "One-quarter of the units", "Half of the units", "All of the units")
mod2$modelName<-"Renters"
names(mod2)<-c("affordability","est","se","uci","lci","Quintile","modelName")
modFrame<-data.frame(rbind(mod1,mod2))
modFrame$Quintile <- factor(modFrame$Quintile, levels = c("None of the units", "One-quarter of the units", "Half of the units", "All of the units"))
interval1<-qnorm((1-.9)/2)
interval2<-qnorm((1-.95)/2)
modFrame$est<-as.numeric(as.character(modFrame$est))
modFrame$se<-as.numeric(as.character(modFrame$se))

modFrame$modelName<-factor(modFrame$modelName, levels=c("Homeowners","Renters"))
modFrame
renters_type_nimby<-ggplot(modFrame, aes(colour=modelName, shape=modelName))+ scale_y_continuous(limits = c(-.2, .2))
renters_type_nimby<-renters_type_nimby+theme_bw()+scale_colour_grey(end=.5)+geom_hline(yintercept=0, colour=gray(1/2), lty=2)
renters_type_nimby<-renters_type_nimby+geom_linerange(aes(x=Quintile, ymin=est-se*interval1, 
                                                          ymax=est+se*interval1), lwd=1, position=position_dodge(width=1/2))
renters_type_nimby<-renters_type_nimby+geom_pointrange(aes(x=Quintile, y=est, ymin=est-se*interval2,
                                                           ymax=est+se*interval2), lwd=1/2,
                                                       position=position_dodge(width=1/2),  fill="WHITE")
renters_type_nimby<-renters_type_nimby+coord_flip()+ labs(y="Change in Probability Building Preferred",x="Affordability Level")
renters_type_nimby<-renters_type_nimby+theme(legend.title=element_blank())+ theme(aspect.ratio = .5)
renters_type_nimby<-renters_type_nimby+theme(plot.margin=unit(c(0,0,0,0),"mm")) + ggtitle("Proximity by Affordability") + guides(colour=guide_legend(reverse=T), shape=guide_legend(reverse=T))
pdf("plots/FigC3_ownership_affordability_effects.pdf", width=6, height=3)
print(renters_type_nimby)
dev.off()

tiff("plots/FigC3_ownership_affordability_effects.tif", res=600, compression="lzw", height=3, width=6, units="in")
print(ideology_affordable)
dev.off()

# Figure C.4 Renters Proximity by each level of affordability, by City ####
quantile(conjoint4$zri_city, probs=seq(0,1,.2), na.rm=T)

zri_values<-c(0,1217,1480,1936,2427,7500)
est1<-rep(NA, length(zri_values))
se1<-rep(NA, length(zri_values))
for(i in 1:5){
  mod1<-amce(select ~ distance+ community + height + site + tenant + units,
             data=subset(renters, zri_city>zri_values[i] & zri_city<=zri_values[i+1]&affordable=="None of the units"), cluster=T, respondent.id = "CaseID")
  est1[i]<-summary(mod1)$amce[5,3]
  se1[i]<-summary(mod1)$amce[5,4]
}
mod1ests<-as.data.frame(cbind(zri_values,est1,se1))
mod1ests$uCI<-est1+se1*1.96
mod1ests$lCI<-est1-se1*1.96

est2<-rep(NA, length(zri_values))
se2<-rep(NA, length(zri_values))
for(i in 1:5){
  mod2<-amce(select ~ distance+ community + height + site + tenant + units,
             data=subset(renters, zri_city>zri_values[i] & zri_city<=zri_values[i+1]&affordable=="One-quarter of the units"), cluster=T, respondent.id = "CaseID")
  est2[i]<-summary(mod2)$amce[5,3]
  se2[i]<-summary(mod2)$amce[5,4]
}
mod2ests<-as.data.frame(cbind(zri_values,est2,se2))
mod2ests$uCI<-est2+se2*1.96
mod2ests$lCI<-est2-se2*1.96

est3<-rep(NA, length(zri_values))
se3<-rep(NA, length(zri_values))
for(i in 1:5){
  mod3<-amce(select ~ distance+ community + height + site + tenant + units,
             data=subset(renters, zri_city>zri_values[i] & zri_city<=zri_values[i+1]&affordable=="Half of the units"), cluster=T, respondent.id = "CaseID")
  est3[i]<-summary(mod3)$amce[5,3]
  se3[i]<-summary(mod3)$amce[5,4]
}
mod3ests<-as.data.frame(cbind(zri_values,est3,se3))
mod3ests$uCI<-est3+se3*1.96
mod3ests$lCI<-est3-se3*1.96

est4<-rep(NA, length(zri_values))
se4<-rep(NA, length(zri_values))
for(i in 1:5){
  mod4<-amce(select ~ distance+ community + height + site + tenant + units,
             data=subset(renters, zri_city>zri_values[i] & zri_city<=zri_values[i+1]&affordable=="All of the units"), cluster=T, respondent.id = "CaseID")
  est4[i]<-summary(mod4)$amce[5,3]
  se4[i]<-summary(mod4)$amce[5,4]
}
mod4ests<-as.data.frame(cbind(zri_values,est4,se4))
mod4ests$uCI<-est4+se4*1.96
mod4ests$lCI<-est4-se4*1.96

#combine data
mod1<-mod1ests[-6,]
mod1$quintle<-c("Least Expensive","2nd","3rd","4th","Most Expensive")
mod1$modelName<-"None of Units"
names(mod1)<-c("zri","est","se","uci","lci","Quintile","modelName")

mod2<-mod2ests[-6,]
mod2$quintle<-c("Least Expensive","2nd","3rd","4th","Most Expensive")
mod2$modelName<-"Quarter of Units"
names(mod2)<-c("zri","est","se","uci","lci","Quintile","modelName")

mod3<-mod3ests[-6,]
mod3$quintle<-c("Least Expensive","2nd","3rd","4th","Most Expensive")
mod3$modelName<-"Half of Units"
names(mod3)<-c("zri","est","se","uci","lci","Quintile","modelName")

mod4<-mod4ests[-6,]
mod4$quintle<-c("Least Expensive","2nd","3rd","4th","Most Expensive")
mod4$modelName<-"All of Units"
names(mod4)<-c("zri","est","se","uci","lci","Quintile","modelName")

modFrame<-data.frame(rbind(mod1,mod2,mod3,mod4))
modFrame$Quintile <- factor(modFrame$Quintile, levels = c("Least Expensive","2nd","3rd","4th","Most Expensive"))
interval1<-qnorm((1-.9)/2)
interval2<-qnorm((1-.95)/2)
modFrame
modFrame$modelName<-factor(modFrame$modelName, levels=c("None of Units","Quarter of Units","Half of Units","All of Units"))
renters_type_nimby<-ggplot(modFrame, aes(colour=modelName, shape=modelName))+ scale_y_continuous(limits = c(-.3, .3))#+ scale_colour_discrete(name="Affordability") +scale_shape_discrete(name="Affordability")
renters_type_nimby<-renters_type_nimby+theme_bw()+theme(legend.title = element_blank())+scale_colour_grey(end=.5)+geom_hline(yintercept=0, colour=gray(1/2), lty=2)
renters_type_nimby<-renters_type_nimby+geom_linerange(aes(x=Quintile, ymin=est-se*interval1, 
                                                          ymax=est+se*interval1), lwd=1, position=position_dodge(width=1/2))
renters_type_nimby<-renters_type_nimby+geom_pointrange(aes(x=Quintile, y=est, ymin=est-se*interval2,
                                                           ymax=est+se*interval2), lwd=1/2,
                                                       position=position_dodge(width=1/2),  fill="WHITE")
renters_type_nimby<-renters_type_nimby+coord_flip()+ labs(y="Change in Probability Building Preferred",x="Average City Rent (Quintiles)")
renters_type_nimby<-renters_type_nimby+theme(plot.margin=unit(c(0,0,0,0),"mm"))+ggtitle(("Renters, Proximity by Affordability and Average Rent by City"))+guides(colour=guide_legend(reverse=T), shape=guide_legend(reverse=T))
pdf("plots/FigC4_renters_type_nimby_quintile_zri_city_all.pdf", width=6, height=3)
print(renters_type_nimby)
dev.off()

tiff("plots/FigC4_renters_type_nimby_quintile_zri_city_all.tif", res=600, compression="lzw", height=3, width=6, units="in")
print(renters_type_nimby)
dev.off()


# Figure C.5 Renters Proximity by Affordability and Average Rent by ZIP ####
quantile(conjoint4$zri, probs=seq(0,1,.2), na.rm=T)
zri_city_values<-c(0,1204,1526,1959,2488,13000)
est1<-rep(NA, length(zri_city_values))
se1<-rep(NA, length(zri_city_values))
for(i in 1:5){
  mod1<-amce(select ~ distance+ community + height + site + tenant + units,
             data=subset(renters, zri>zri_city_values[i] & zri<=zri_city_values[i+1]&luxury==1), cluster=T, respondent.id = "CaseID")
  est1[i]<-summary(mod1)$amce[5,3]
  se1[i]<-summary(mod1)$amce[5,4]
}
mod1ests<-as.data.frame(cbind(zri_city_values,est1,se1))
mod1ests$uCI<-est1+se1*1.96
mod1ests$lCI<-est1-se1*1.96
est2<-rep(NA, length(zri_city_values))
se2<-rep(NA, length(zri_city_values))
for(i in 1:5){
  mod2<-amce(select ~ distance+ community + height + site + tenant + units,
             data=subset(renters, zri>zri_city_values[i] & zri<=zri_city_values[i+1]&luxury==0), cluster=T, respondent.id = "CaseID")
  est2[i]<-summary(mod2)$amce[5,3]
  se2[i]<-summary(mod2)$amce[5,4]
}
mod2ests<-as.data.frame(cbind(zri_city_values,est2,se2))
mod2ests$uCI<-est2+se2*1.96
mod2ests$lCI<-est2-se2*1.96
#combine data
mod1<-mod1ests[-6,]
mod1$quintle<-c("Least Expensive","2nd","3rd","4th","Most Expensive")
mod1$modelName<-"Market Rate"
names(mod1)<-c("zri","est","se","uci","lci","Quintile","modelName")
mod2<-mod2ests[-6,]
mod2$quintle<-c("Least Expensive","2nd","3rd","4th","Most Expensive")
mod2$modelName<-"Affordable"
names(mod2)<-c("zri","est","se","uci","lci","Quintile","modelName")
modFrame<-data.frame(rbind(mod1,mod2))
modFrame$Quintile <- factor(modFrame$Quintile, levels = c("Least Expensive", "2nd","3rd", "4th", "Most Expensive"))
interval1<-qnorm((1-.9)/2)
interval2<-qnorm((1-.95)/2)
modFrame
modFrame$modelName<-factor(modFrame$modelName, levels=c("Market Rate","Affordable"))
renters_type_nimby<-ggplot(modFrame, aes(colour=modelName, shape=modelName))+ scale_y_continuous(limits = c(-.25, .25))
renters_type_nimby<-renters_type_nimby+theme_bw()+scale_colour_grey(end=.5)+geom_hline(yintercept=0, colour=gray(1/2), lty=2)
renters_type_nimby<-renters_type_nimby+geom_linerange(aes(x=Quintile, ymin=est-se*interval1, 
                                                          ymax=est+se*interval1), lwd=1, position=position_dodge(width=1/2))
renters_type_nimby<-renters_type_nimby+geom_pointrange(aes(x=Quintile, y=est, ymin=est-se*interval2,
                                                           ymax=est+se*interval2), lwd=1/2,
                                                       position=position_dodge(width=1/2),  fill="WHITE")
renters_type_nimby<-renters_type_nimby+coord_flip()+ labs(y="Change in Probability Building Preferred",x="Average ZIP Rent (Quintiles)")
renters_type_nimby<-renters_type_nimby+theme(legend.title=element_blank()) + theme(aspect.ratio = .5)
renters_type_nimby<-renters_type_nimby+theme(plot.margin=unit(c(0,0,0,0),"mm"))+ ggtitle("Renters, Proximity by Affordability and Average Rent by ZIP") +guides(colour=guide_legend(reverse=T), shape=guide_legend(reverse=T))
pdf("plots/FigC5_renters_type_nimby_ZIP.pdf", width=6, height=3)
print(renters_type_nimby)
dev.off()

tiff("plots/FigC5_renters_type_nimby_ZIP.tif", res=600, compression="lzw", height=3, width=6, units="in")
print(renters_type_nimby)
dev.off()


# Figure C.6 Homeowners Proximity by each level of affordability, by City ####
quantile(conjoint4$zri_city, probs=seq(0,1,.2), na.rm=T)
zri_values<-c(0,1217,1480,1936,2427,7500)
est1<-rep(NA, length(zri_values))
se1<-rep(NA, length(zri_values))
for(i in 1:5){
  mod1<-amce(select ~ distance+ community + height + site + tenant + units,
             data=subset(owners, zri_city>zri_values[i] & zri_city<zri_values[i+1]&affordable=="None of the units"), cluster=T, respondent.id = "CaseID")
  est1[i]<-summary(mod1)$amce[5,3]
  se1[i]<-summary(mod1)$amce[5,4]
}
mod1ests<-as.data.frame(cbind(zri_values,est1,se1))
mod1ests$uCI<-est1+se1*1.96
mod1ests$lCI<-est1-se1*1.96


est2<-rep(NA, length(zri_values))
se2<-rep(NA, length(zri_values))
for(i in 1:5){
  mod2<-amce(select ~ distance+ community + height + site + tenant + units,
             data=subset(owners, zri_city>zri_values[i] & zri_city<zri_values[i+1]&affordable=="One-quarter of the units"), cluster=T, respondent.id = "CaseID")
  est2[i]<-summary(mod2)$amce[5,3]
  se2[i]<-summary(mod2)$amce[5,4]
}
mod2ests<-as.data.frame(cbind(zri_values,est2,se2))
mod2ests$uCI<-est2+se2*1.96
mod2ests$lCI<-est2-se2*1.96

est3<-rep(NA, length(zri_values))
se3<-rep(NA, length(zri_values))
for(i in 1:5){
  mod3<-amce(select ~ distance+ community + height + site + tenant + units,
             data=subset(owners, zri_city>zri_values[i] & zri_city<zri_values[i+1]&affordable=="Half of the units"), cluster=T, respondent.id = "CaseID")
  est3[i]<-summary(mod3)$amce[5,3]
  se3[i]<-summary(mod3)$amce[5,4]
}
mod3ests<-as.data.frame(cbind(zri_values,est3,se3))
mod3ests$uCI<-est3+se3*1.96
mod3ests$lCI<-est3-se3*1.96

est4<-rep(NA, length(zri_values))
se4<-rep(NA, length(zri_values))
for(i in 1:5){
  mod4<-amce(select ~ distance+ community + height + site + tenant + units,
             data=subset(owners, zri_city>zri_values[i] & zri_city<zri_values[i+1]&affordable=="All of the units"), cluster=T, respondent.id = "CaseID")
  est4[i]<-summary(mod4)$amce[5,3]
  se4[i]<-summary(mod4)$amce[5,4]
}
mod4ests<-as.data.frame(cbind(zri_values,est4,se4))
mod4ests$uCI<-est4+se4*1.96
mod4ests$lCI<-est4-se4*1.96

#combine data
mod1<-mod1ests[-6,]
mod1$quintle<-c("Least Expensive","2nd","3rd","4th","Most Expensive")
mod1$modelName<-"None of Units"
names(mod1)<-c("zri","est","se","uci","lci","Quintile","modelName")

mod2<-mod2ests[-6,]
mod2$quintle<-c("Least Expensive","2nd","3rd","4th","Most Expensive")
mod2$modelName<-"Quarter of Units"
names(mod2)<-c("zri","est","se","uci","lci","Quintile","modelName")

mod3<-mod3ests[-6,]
mod3$quintle<-c("Least Expensive","2nd","3rd","4th","Most Expensive")
mod3$modelName<-"Half of Units"
names(mod3)<-c("zri","est","se","uci","lci","Quintile","modelName")

mod4<-mod4ests[-6,]
mod4$quintle<-c("Least Expensive","2nd","3rd","4th","Most Expensive")
mod4$modelName<-"All of Units"
names(mod4)<-c("zri","est","se","uci","lci","Quintile","modelName")

modFrame<-data.frame(rbind(mod1,mod2,mod3,mod4))
modFrame$Quintile <- factor(modFrame$Quintile, levels = c("Least Expensive","2nd","3rd","4th","Most Expensive"))
interval1<-qnorm((1-.9)/2)
interval2<-qnorm((1-.95)/2)
modFrame
modFrame$modelName<-factor(modFrame$modelName, levels=c("None of Units","Quarter of Units","Half of Units","All of Units"))
renters_type_nimby<-ggplot(modFrame, aes(colour=modelName, shape=modelName))+ scale_y_continuous(limits = c(-.3, .3))#+ scale_colour_discrete(name="Affordability") +scale_shape_discrete(name="Affordability")
renters_type_nimby<-renters_type_nimby+theme_bw()+theme(legend.title = element_blank())+scale_colour_grey(end=.5)+geom_hline(yintercept=0, colour=gray(1/2), lty=2)
renters_type_nimby<-renters_type_nimby+geom_linerange(aes( x=Quintile, ymin=est-se*interval1, 
                                                           ymax=est+se*interval1), lwd=1, position=position_dodge(width=1/2))
renters_type_nimby<-renters_type_nimby+geom_pointrange(aes(x=Quintile, y=est, ymin=est-se*interval2,
                                                           ymax=est+se*interval2), lwd=1/2,
                                                       position=position_dodge(width=1/2),  fill="WHITE")
renters_type_nimby<-renters_type_nimby+coord_flip()+ labs(y="Change in Probability Building Preferred", x="Average Rent by City (Quintiles)")
renters_type_nimby<-renters_type_nimby+theme(plot.margin=unit(c(0,0,0,0),"mm"))+ggtitle(("Homeowners, Proximity by Affordability and Average Rent by City"))+guides(colour=guide_legend(reverse=T), shape=guide_legend(reverse=T))
pdf("plots/FigC6_owners_type_nimby_quintile_zri_city.pdf", width=6, height=3)
print(renters_type_nimby)
dev.off()

tiff("plots/FigC6_owners_type_nimby_quintile_zri_city.tif", res=600, compression="lzw", height=3, width=6, units="in")
print(renters_type_nimby)
dev.off()


######################
#### SF EXIT POLL 
######################

setwd("/Users/Michael/Dropbox/apsr")

#bring in data
final<-read.csv("data/sfDataAPSR.csv", stringsAsFactors = F)

owners<-subset(final, ownership_dummy==1)
renters<-subset(final, ownership_dummy==0)

#experiment randomization
control<-subset(final, version==5 | version==6)
control_owners<-subset(control, ownership_dummy==1)
control_renters<-subset(control, ownership_dummy==0)

#percent support
table(control_owners$ten_plan_dummy)
87/(87+33) #.73
table(control_renters$ten_plan_dummy)
171/(32+171) #.84


# Figure 1. Support for Micro-scale Ban by Support for Macro-scale Supply ####
control_owners_yes<-subset(control_owners, ten_plan_dummy==1)
control_owners_no<-subset(control_owners, ten_plan_dummy==0)

control_renters_yes<-subset(control_renters, ten_plan_dummy==1)
control_renters_no<-subset(control_renters, ten_plan_dummy==0)

table(control_owners_yes$prop_i_ban_dummy)
29/(29+49) #.37
table(control_owners_no$prop_i_ban_dummy)
15/30 #.50

table(control_renters_yes$prop_i_ban_dummy)
79/(79+73) #.52
table(control_renters_no$prop_i_ban_dummy)
22/(22+5) #.81

support<-c(.37, .50, .52, .81)
Tenant<-c("Homeowners", "Homeowners", "Renters", "Renters")
supply<-c("Pro-Supply", "Anti-Supply","Pro-Supply","Anti-Supply")
supply<-factor(supply, levels=c("Pro-Supply", "Anti-Supply"))
ban_plot<-data.frame(Tenant, support, supply)

# Figure
exitpoll_ban<-ggplot(data=ban_plot, aes(x=Tenant, y=support, fill=Tenant))+
  geom_bar(colour="black",stat="identity", position=position_dodge()) + facet_wrap(~supply) +
  ylab("Support (Percent)") + theme(legend.position="none") +ggtitle( "Support for Micro-scale Ban by Support for Macro-scale Supply")+
  theme_bw()+scale_fill_grey()
exitpoll_ban
pdf("plots/Fig1_exitpoll_ban.pdf", height=3, width=6)
plot(exitpoll_ban)
dev.off()

tiff("plots/Fig1_exitpoll_ban.tif", res=600, compression="lzw", height=3, width=6, units="in")
print(exitpoll_ban)
dev.off()


# Table A.3: Policy Proposals, SF ####
#model control supply support
simple_control<-(lm(ten_plan_dummy ~ ownership_dummy, final)); summary(simple_control)
simple_control_se<-sqrt(diag(vcovHC(simple_control, type="HC1")))

full_control<-(lm(ten_plan_dummy ~ ownership_dummy  +  scale(ideology_num) +scale(income_num) + white_dummy  + age+ gender_dummy , control)); summary(full_control)
full_control_se<-sqrt(diag(vcovHC(full_control, type="HC1")))

# Table - Control Condition
#Supplementary Data
stargazer(simple_control, full_control, title="Ten Percent Supply Increase, San Francisco",  label="ten_pct_dummy",
          dep.var.labels=c("Support Supply Increase"), dep.var.labels.include = F, dev.var.caption="",
          column.labels=c("Bivariate", "Full"),
          covariate.labels=c("Homeownership","Ideology","Income, Log","White, Non-Hispanic","Age","Male"),
          omit.stat=c("ser","f"), digits=2, align=T,
          initial.zero = F, font.size="small", star.cutoffs = NA, omit.table.layout="n",
          se=list(simple_control_se, full_control_se), no.space=T, omit=c("name"))

#Ban within all
table(owners$prop_i_ban_dummy)
187/(281+187) #.40

table(renters$prop_i_ban_dummy)
508/(318+508) #.62

#model ban
simple_prop_i_ban<-(lm(prop_i_ban_dummy ~ ownership_dummy  , final)); summary(simple_prop_i_ban)
simple_prop_i_ban_se<-sqrt(diag(vcovHC(simple_prop_i_ban, type="HC1")))

full_prop_i_ban<-(lm(prop_i_ban_dummy ~ ownership_dummy  +  scale(ideology_num) +scale(income_num) + white_dummy  + age+ gender_dummy , final)); summary(full_prop_i_ban)
full_prop_i_ban_se<-sqrt(diag(vcovHC(full_prop_i_ban, type="HC1")))

#Table. Ban support
stargazer(simple_prop_i_ban, full_prop_i_ban, title="Neighborhood Ban, San Francisco",  label="prop_i_ban_dummy",
          dep.var.labels=c("Support Supply Increase"), dep.var.labels.include = F, dev.var.caption="",
          column.labels=c("Bivariate", "Full"),
          covariate.labels=c("Homeownership","Ideology","Income, Log","White, Non-Hispanic","Age","Male"),
          omit.stat=c("ser","f"), digits=2, align=T,
          initial.zero = F, font.size="small", star.cutoffs = NA, omit.table.layout="n",
          se=list(simple_prop_i_ban_se, full_prop_i_ban_se), no.space=T, omit=c("name"))

# Table. Combine these two tables
stargazer(simple_control, full_control, simple_prop_i_ban, full_prop_i_ban, title="Policy Proposals, San Francisco Sample",  label="sf_policies",
          dep.var.labels.include = F, dev.var.caption="",
          column.labels=c("10 Pct Supply", "NIMBY Ban Proposal" ), column.separate = c(2, 2),
          covariate.labels=c("Homeownership","Ideology","Income, Log","White, Non-Hispanic","Age","Male"),
          omit.stat=c("ser","f"), digits=2, align=T, type="text",
          initial.zero = F, font.size="small", star.cutoffs = NA, omit.table.layout="n",
          se=list(simple_control_se, full_control_se, simple_prop_i_ban_se, full_prop_i_ban_se), no.space=T )


########################
#### Recontacted SF
########################

setwd("/Users/Michael/Dropbox/apsr")
socpoc_sf<-read.csv("data/socpocSFAPSR.csv")

#assign groups
owners<-subset(socpoc_sf, own==1 )
renters<-subset(socpoc_sf, own==0)


#### Rent Control Tests 
#proposition a
summary(lm(prop_a_dummy ~ rent_control + income_cat + ideology+ whitenh, renters))

#proposition d
summary(lm(prop_d_dummy ~ rent_control + income_cat + ideology+ whitenh, renters))

#proposition f
summary(lm(prop_f_dummy ~ rent_control + income_cat + ideology+ whitenh, renters))

#proposition i
summary(lm(prop_i_dummy ~ rent_control + income_cat + ideology+ whitenh, renters))

#neighborhood ban
summary(lm(prop_i_ban_dummy ~ rent_control + income_cat + ideology+ whitenh, renters))

#housing
summary(lm(supply_dummy ~ rent_control + income_cat + ideology+ whitenh, renters))


#########################################
##### Recontacted SF - Conjoint Data 
#########################################

conjoint_sf<-read.csv("data/conjointSFAPSR.csv")

conjoint_sf$distance <- factor(conjoint_sf$distance,levels= c("2 miles (40 minute walk)", "1 mile (20 minute walk)", "1/2 mile (10 minute walk)",
                                                              "1/8 mile (2 minute walk)"))
conjoint_sf$community <- factor(conjoint_sf$community,levels= c("No opinion", "Support the building", "Oppose the building"))
conjoint_sf$affordable <- factor(conjoint_sf$affordable,levels= c("None of the units", "One-quarter of the units", "Half of the units", "All of the units"))
conjoint_sf$height <- factor(conjoint_sf$height,levels= c("2 stories", "3 stories", "6 stories", "12 stories"))
conjoint_sf$site <- factor(conjoint_sf$site, levels=c("Empty building","Parking lot","Open field","Historically-designated building"))

#add indicators 
conjoint_sf$city_interest_low<-as.factor(ifelse(conjoint_sf$city_interest<0,1,0))
conjoint_sf$prop_i_ban_dummy<-as.factor(ifelse(conjoint_sf$prop_i_ban_dummy==1,1,0))
conjoint_sf$luxury<-as.factor(ifelse(conjoint_sf$affordable=="None of the units",1,0))

#subgroups
owners<-subset(conjoint_sf, own==1)
renters<-subset(conjoint_sf, own==0)

owners_aff<-subset(owners, luxury==0)
owners_lux<-subset(owners, luxury==1)

renters_aff<-subset(renters, luxury==0)
renters_lux<-subset(renters, luxury==1)

# Figure A.1. Recontacted Conjoint San Francisco Sample ####
renters_prop_i_yes_mod<-amce(selected ~ distance  + community + height + site + tenant + units,
                             data=subset(renters_lux, prop_i_ban_dummy==1&supply_dummy==1),cluster=T,  respondent.id = "respondent")
renters_prop_i_no_mod<-amce(selected ~ distance  + community + height + site + tenant + units,
                            data=subset(renters_lux, prop_i_ban_dummy==0&supply_dummy==1), cluster=T,  respondent.id = "respondent")
renters_prop_i_yes_mod_frame<-data.frame(Variable=(summary(renters_prop_i_yes_mod)$amce)$Level,
                                         Coefficient = (summary(renters_prop_i_yes_mod)$amce)$Estimate,
                                         SE=(summary(renters_prop_i_yes_mod)$amce)$'Std. Err',
                                         modelName="Supporters")
renters_prop_i_no_mod_frame<-data.frame(Variable=(summary(renters_prop_i_no_mod)$amce)$Level,
                                        Coefficient = (summary(renters_prop_i_no_mod)$amce)$Estimate,
                                        SE=(summary(renters_prop_i_no_mod)$amce)$'Std. Err',
                                        modelName="Opponents")
rentersPropIFrame<-data.frame(rbind(renters_prop_i_yes_mod_frame, renters_prop_i_no_mod_frame))
rentersPropIFrame<-subset(rentersPropIFrame, Variable=="1/8 mile (2 minute walk)"|Variable=="1/2 mile (10 minute walk)"|Variable=="1 mile (20 minute walk)")
rentersPropIFrameIntercepts<-data.frame(Variable=c("2 miles (40 minute walk)", "2 miles (40 minute walk)") ,Coefficient=c(0,0), SE=c(0,0), modelName=c("Supporters", "Opponents"))
rentersPropIFrame<-data.frame(rbind(rentersPropIFrame,rentersPropIFrameIntercepts))
interval1<-qnorm((1-.9)/2)
interval2<-qnorm((1-.95)/2)
rentersPropIFrame$Variable <- factor(rentersPropIFrame$Variable, levels = c("1/8 mile (2 minute walk)","1/2 mile (10 minute walk)", "1 mile (20 minute walk)", "2 miles (40 minute walk)"))
rentersPropIFrame
rentersPropINimby<-ggplot(rentersPropIFrame, aes(colour=modelName, shape=modelName))+scale_shape_manual(values=c(4,5))
#+  scale_color_manual(values=c('#FF9933','#3399FF'))
rentersPropINimby<-rentersPropINimby+theme_bw()+scale_colour_grey(end=.5)+geom_hline(yintercept=0, colour=gray(1/2), lty=2)
rentersPropINimby<-rentersPropINimby+geom_linerange(aes(x=Variable, ymin=Coefficient-SE*interval1, 
                                                        ymax=Coefficient+SE*interval1), lwd=1, position=position_dodge(width=1/2))
rentersPropINimby<-rentersPropINimby+geom_pointrange(aes(x=Variable, y=Coefficient, ymin=Coefficient-SE*interval2,
                                                         ymax=Coefficient+SE*interval2), lwd=1/2,
                                                     position=position_dodge(width=1/2), fill="WHITE")
rentersPropINimby<-rentersPropINimby+coord_flip()+labs(y="Change in Probability Building Preferred")
rentersPropINimby<-rentersPropINimby+theme(legend.title=element_blank(),axis.title.y=element_blank())+ theme(aspect.ratio = .5)
rentersPropINimby<-rentersPropINimby+theme(plot.margin=unit(c(0,0,0,0),"mm")) + ggtitle("Renters, Proximity and Ban Support (Market Rate)") + guides(colour=guide_legend(reverse=T), shape=guide_legend(reverse=T))

pdf("plots/FigA1_renters_prop_i_ban_nimby.pdf", width=6, height=3)
print(rentersPropINimby)
dev.off()

tiff("plots/FigA1_renters_prop_i_ban_nimby.tif", res=600, compression="lzw", height=3, width=6, units="in")
print(rentersPropINimby)
dev.off()


###################################################
#### Descriptive Statistics for Sampling Frame ####
###################################################

#ZIP Code Characteristics
zips<-read.csv("data/hankinsonzips.csv")
zips$ZIP<-formatC(zips$ZIP, width=5, format="d", flag="0")

#pull in all census data for zips? then subset?
#percent white
race<-read.csv("data/acs2015race/ACS_15_5YR_DP05_with_ann.csv", stringsAsFactors = F)
race1<-race[,c(1:4,8,290)]
colnames(race1)<-c("id", "id2", "zcta","pop","male","whitenh")
race2<-race1[-1,]

#income, household 2015 ACS 5-year
income<-read.csv("data/acs2015income/ACS_15_5YR_S1902_with_ann.csv", stringsAsFactors = F)
income1<-income[,c(1:3,8)]
colnames(income1)<-c("id", "id2", "zcta","income_mean")
income2<-income1[-1,]

incomeMedian<-read.csv("data/acs15incomemedian/ACS_15_5YR_S1903_with_ann.csv", stringsAsFactors = F)
incomeMedian1<-incomeMedian[,c(1:3,6)]
colnames(incomeMedian1)<-c("id", "id2", "zcta","income_median")
incomeMedian2<-incomeMedian1[-1,]

#homeownership, percent 2015 ACS 5-year
housing<-read.csv("data/acs2015housing/ACS_15_5YR_DP04_with_ann.csv", stringsAsFactors = F)
housing1<-housing[,c(1:3,186)]
colnames(housing1)<-c("id", "id2", "zcta","owner_occ")
housing2<-housing1[-1,]

#merge acs estimates
acs_zips<-merge(race2, incomeMedian2, by="id")
acs_zips1<-merge(acs_zips, housing2, by="id")

#subset to my 4068.
acs_zips2<-merge(zips, acs_zips1, by.x="ZIP", by.y="id2")
acs_zips2<-subset(acs_zips2,pop>0)

sum(as.numeric(acs_zips2$whitenh)*as.numeric(acs_zips2$pop))/(sum(as.numeric(acs_zips2$pop)))
# .463
sum(as.numeric(acs_zips2$owner_occ)*as.numeric(acs_zips2$pop), na.rm=T)/(sum(as.numeric(acs_zips2$pop)))
# .495
sum(as.numeric(acs_zips2$income_median)*as.numeric(acs_zips2$pop), na.rm=T)/(sum(as.numeric(acs_zips2$pop)))
# 57,101
sum(as.numeric(acs_zips2$male))/sum(as.numeric(acs_zips2$pop))
#.489





